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Abstract. A kinetic model of particle acceleration in su- 
pernova remnants (SNRs) is extended to study the cos- 
mic ray (CR) and associated high-energy 7-ray production 
during SN shock propagation through the inhomogeneous 
circumstellar medium of a progenitor star that emits a 
wind. The wind forms a low-density bubble surrounded by 
a swept-up shell of interstellar matter. 7-rays are produced 
as a result of decay of pions which in turn are the result of 
collisions of CRs with nuclei of the thermal plasma. The 
time evolution of the SNRs is followed numerically, taking 
into account the nonlinear backreaction of the accelerated 
CRs. The model for SNRs includes injection of suprather- 
mal particles at the shock front and heating of the thermal 
plasma due to dissipation of Alfven waves in the precursor 
region. Examples typical for SN type lb and SN type II 
explosions are considered. Apart from the confirmation of 
the known fact that acceleration is extremely rapid and 
that the upper momentum cutoff is reached almost imme- 
diately after the explosion due to the high wind magnetic 
field, it is also shown that the CRs are accelerated with 
a high efficiency. Depending on the circumstellar parame- 
ters, 20% to 40% of the SN explosion energy is absorbed 
by CRs during the SNR evolution for a moderate injec- 
tion rate, when a fraction 77 = 10"'^ of the gas particles, 
swept up by the supernova shock, is accelerated. The CR 
momentum spectrum, ultimately produced in the SNRs, 
has a power law form N oc p^"' with an index 7 — 2.0 
to 2.1 in a wide energy range up to a maximum energy, 
which is about 10^^ eV, if the CR diffusion coefficient is 
as small as the Bohm limiting value. It is to be expected 
that the resulting CR chemical composition at high ener- 
gies reflects more the stellar wind composition, whereas at 
lower energies it corresponds more to that of the average 
interstellar medium. The expected Tr'^-decay 7-ray flux is 
however considerably lower than in the case of a uniform 
interstellar medium; a relatively high 7-ray luminosity in 
the band e-y ^ 1 TeV, detectable at distances of several kpc, 
is only expected in the case of a relatively dense ISM with 



a number density above 10 cm~'^. Extremely high 7-ray 
emission may be produced when the SN shock propagates 
through the slow dense wind of a red supergiant, the pro- 
genitor of a SN type II. In this case SNRs might be visible 
in 7-rays for several hundred years out to distances of tens 
of kpc. For the case of a SN type lb the expected Tr'^-decay 
7-ray TeV-energy flux during the whole SNR evolution 



ber density is less than 0.3 cm ^ . 
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1. Introduction 

Considerable efforts have been made during the last years 
to empirically confirm the theoretical expectation that the 
main part of the Galactic cosmic rays (CRs) originates in 
supernova remnants (SNRs) . Theoretically progress in the 
solution of this problem has been due to the development 
of the theory of diffusive shock acceleration (see, for ex- 
ample, reviews by Drury 1983; Blandford & Eichler 1987; 
Berezhko & Krymsky 1988). Although still incomplete, 
the theory is able to explain the main characteristics of the 
observed CR spectrum under several reasonable assump- 
tions, at least up to an energy of 10"'^^ 10"'^^ eV. Direct 
information about the dominant nucleonic CR component 
in SNRs can only be vobtained from 7-ray observations. If 
this nuclear component is strongly enhanced inside SNRs 
then through inelastic nuclear collisions, leading to pion 
production and subsequent decay, 7-rays will be produced. 

CR acceleration in SNRs expanding in a uniform in- 
terstellar medium (ISM) (Drury et al. 1989; Markiewicz et 
al. 1990; Dorfi 1990), and the properties of the associated 
7-ray emission (Dorfi 1991; Drury et al. 1994) were investi- 
gated in a number of studies (we mention here only those 
papers which include the effects of shock geometry and 
time-dependent nonlinear CR backreaction; for a review 
of others which deal with the test particle approximation. 
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sec for example Drury 1983; Blandford & Eichler 1987; 
Berezhko & Krymsky 1988; Volk 1997). All of these stud- 
ies are based on a two-fluid hydro dynamical approach and 
directly employ the assumption that the expanding SN 
shock is locally plane; as dynamic variables for the CRs 
the pressure and the energy density are determined. Their 
characteristics are sometimes essentially different from the 
results obtained in a kinetic approach (Kang & Jones 1991; 
Berezhko et al. 1994, 1995, 1996) which consistently takes 
the role of shock geometry and nonlinear CR backreaction 
into account. First of all, in kinetic theory the form of the 
spectrum of accelerated CRs and their maximum energy 
are calculated selfconsistently. In particular, the maximum 
particle energy Cmax, achieved at any given evolution- 
ary stage, is determined by geometrical factors (Berezhko 
1996), in contrast to the hydrodynamic models which in 
fact postulate that the value of emax{t) is determined by 
the time interval t that has passed since the explosion 
(Drury et al. 1989; Markiewicz et al. 1990; Dorfi 1990; 
Jones & Kang 1992). Although the difference between the 
values of Cmax in the two cases is not very large, it criti- 
cally influences the structure and evolution of the shock. 
For example, the shock never becomes completely modi- 
fled (smoothed) by the CR backreaction (Kang & Jones 
1991; Berezhko ct al. 1994, 1995, 1996). Together with the 
smooth precursor, the shock transition always contains a 
relatively strong subshock which heats the swept-up gas 
and leads to the injection of suprathermal gas particles 
into the acceleration process. In this sense diffusive shock 
acceleration is somewhat less efficient than predicted by 
hydrodynamic models. Acceleration always requires some 
freshly injected particles which are generated during gas 
heating. This prediction is in agreement with the observa- 
tions that show significant gas heating in young SNRs. 

On the other hand, the shock modification by the 
CR backreaction is even greater than predicted by hy- 
drodynamic models. The total shock compression ratio 
(7 (X M3/4 (Berezhko et al. 1996) is a monotonically ris- 
ing function of the shock Mach number M and can sig- 
nificantly exceed the classical value 4, or even the value 
a = 7 that corresponds to a postshock medium dominated 
by relativistic CRs. This result is a direct consequence 
of the fact that CR acceleration at the front of a spher- 
ically expanding shock is accompanied by a temporally 
increasing dilution of shock energy in the form of CRs in 
the precursor region (Drury et al. 1995; Berezhko 1996). 
Qualitatively this leads to the same effect as if energetic 
particles or thermal gas energy were radiated away from 
the shock (Berezhko 1996; Berezhko & Elfison 1999). We 
note that this energy dilution effect is much stronger than 
the effect of a decreasing specific heat ratio due to the 
conversion of shock energy into relativistic CRs. 

The additional gas heating in the precursor region due 
to Alfven wave dissipation significantly restricts the shock 
compression ratio (7{t) even though its value still consid- 
erably exceeds 4 for a strong shock (Berezhko et al. 1996; 



Berezhko & Volk 1997). The precursor gas heating has a 
much less pronounced effect in the hydrodynamic descrip- 
tion (e.g. Dorfi 1990). 

As far as the expected 7-ray emission, produced 
in SNRs by the nuclear CR component, is concerned, 
there are less significant differences between the kinetic 
(Berezhko & Volk 1997) and the hydrodynamic (Dorfi 
1991; Drury et al. 1994) predictions, even though these 
differences are not unimportant. The reason is that the 
peak value of the 7-ray fiux is mainly determined by the 
fraction of overall hydrodynamic explosion energy that is 
transferred to CRs, i.e. by the overall efficiency of CR ac- 
celeration which is not strongly dependent on the model 
used. 

The differences in the time variation of the predicted 

7-ray fluxes arc more essential. Kinetic theory (Berezhko 
& Volk 1997) revealed much more effective CR and 7- 
ray production during the free expansion phase. It also 
shows a more rapid decrease of the 7-ray flux during the 
subsequent Sedov phase after reaching its peak value, due 
to the different spatial distributions of thermal gas and 
CRs inside SNRs. This energy-dependent lack of overlap 
is not taken into account in hydrodynamic models. 

The model developed by Ellison and co-workers to de- 
scribe diffusive shock acceleration of CRs and to predict 
the expected 7-ray emission from SNRs, uses a Monte 
Carlo (MC) simulation of pitch angle scattering. It is ap- 
plied in the same way to all particles, both CR and thermal 
particles, and is correspondingly solved numerically from 
the outset (e.g. Ellison et al. 1995; Baring et al. 1999). 
The energetic particle population (i.e. accelerated CRs) 
naturally arises in this approach as a high energy tail of 
the distribution of gas particles which undergo heating 
in the strongest gradient of the selfconsistently modified 
shock transition, the subshock in our notation. In this way 
the model incorporates not only selfconsistent CR accel- 
eration but also gas heating and particle injection into 
the acceleration process. Since it is a plane-wave steady 
state model it includes as a free parameter the value of 
the CR cutoff energy (max, which can not be calculated 
in this kind of model. On the other hand, at sufficiently 
high particle energies the description corresponds to a dif- 
fusive approach based on the diffusive transport equation. 
Therefore our time-dependent kinetic description in spher- 
ical symmetry, with an injection rate that corresponds to 
the MC-model, and the MC-model, in turn incorporat- 
ing the same value of CR cutoff energy, give identical CR 
spectra at any stage of the free expansion phase for the 
same set of ISM and shock parameters (Berezhko & El- 
lison 2000). This means that in its free expansion phase 
the SNR evolution can be represented as a sequence of 
quasi- stationary states; each of them can be reproduced 
in the framework of the plane-wave steady state descrip- 
tion with the corresponding value of the CR cutoff en- 
ergy £max{t)- The situation becomes more complicated in 
the subsequent Sedov phase due to the existence of a so- 
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called escaping particle population near the shock front 
(Berezhko & Ellison 1999b), which is a purely nonstation- 
ary phenomenon (Berezhko 1986a; Berezhko & Krymsky 
1988) and therefore can not be reproduced in the frame- 
work of a steady state approach. As far as the overall CR 
spectrum and the associated 7-ray flux at any given evo- 
lutionary phase are concerned, the plane wave description 
can only give an approximate estimate because it does not 
include such important physical factors as CR adiabatic 
cooling and the incomplete overlap of the spatial distribu- 
tions of gas and CRs inside the expanding SNR. 

The evolution of SNRs in a uniform ISM is typical only 
for SNc type la. SNc of type lb and II, which are more 
numerous in our Galaxy, explode into an inhomogeneous 
circumstellar medium, formed by the intensive wind of 
their massive progenitor stars (Losinskaya 1991). This is 
at least true for progenitor masses M > 20Mq whose main 
sequence winds are very energetic. Their time-integrated 
hydrodynamic energy output is so large that they create 
wind cavities which are about as large as the size of a rem- 
nant of the subsequent Supernova explosion, if it would 
occur in a uniform ISM of the same density. The more fre- 
quent lower mass core collapse progenitors {8<M120Mq) 
have at least a massive Red Supergiant wind region be- 
fore the transition to the external ISM. Therefore the de- 
scription of SNR evolution and CR acceleration becomes 
considerably more diflicult for those cases. Due to the com- 
plicated structure of the medium one can expect that the 
difference between kinetic and hydrodynamic model pre- 
dictions will be even more pronounced than in the case of 
a uniform ISM. 

Up to now, except for the paper by Berezhko & Volk 
(1995) where preliminary results of the present work were 
presented, the consideration of CR acceleration in the case 
of a SN type lb and SN type II was considered only for SN 
shock propagation through the region of the supersonic 
stellar wind (Berezinskii & Ptuskin 1988; Volk & Bier- 
mann 1988; Jones & Kang 1992; Kirk et al. 1995). Since 
this region is characterized by an initially very strong if 
spatially decreasing magnetic field of stellar origin and a 
very high shock speed, the maximum energy of CRs ac- 
celerated in the wind region is achieved extremely rapidly 
during a period of time comparable with the SNR age at 
each phase, and then it remains nearly constant due to 
a balance of acceleration gains and adiabatic losses (Volk 
& Biermann 1988). As it is shown below, the expected 
CR energy spectrum is close to being selfsimilar where 
the amplitude monotonically decreases with time but the 
shape remains invariant (see Appendix). Thus the situa- 
tion is very different from that for a uniform circumstellar 
medium with a standard magnitude ISM field strength, 
where the maximiim CR energy is only achieved at the 
end of the sweep-up phase, after hundreds or even thou- 
sands of years. 

The direct dependence on the stellar surface field also 
raised hopes that the CR acceleration in the wind would 



constitute a major source of the observed Galactic CRs 
with energies e>10^'* eV. The most optimistic assump- 
tions predict CR generation in these regions up to an en- 
ergy of 3 X lO'^^ cV (Biermann 1993). In this context we 
would like to make the following remarks: 

First of all, the supersonic wind region contains a rel- 
atively small amount of mass M < IMq which is typically 
smaller than the ejected mass M^j ~ IOMq. Therefore the 
fraction of the SN explosion energy which can be trans- 
fered into CRs is typically lower than required for the ob- 
served CR spectrum, especially in the case of a SN type 
lb (see below). 

Second, to achieve a particle energy significantly 
greater than 10^^ eV one needs to assume an unusually 
high magnetic field in the wind. This tends to violate the 
condition > Ca- necessary for the existence of a su- 
personic flow in the flrst place, as pointed out by Axford 
(1994), where Vw is the wind speed and Ca is the Alfven 
velocity. 

Third, the wind magnetic fleld is largely azimuthal, 
except in the very polar region near the stellar rotation 
axis, and the outer SNR shock is almost perpendicular. 
Therefore not only the suprathermal ion injection is much 
less efficient but also the CR acceleration itself appears 
not so efficient as in the case of a quasi-parallel shock. Al- 
though one can give some physical arguments which assert 
that particle injection with their subsequent acceleration 
takes place also in this case, one needs a more rigorous 
treatment of this complex question. An example for the 
acceleration at a shock propagating through a stationary 
stochastic large yyscalc magnetic field that is on average 
perpendicular to the shock normal is given by Kirk et al. 
(1996) based on the process of anomalous diffusion (see 
Chuvilgin & Ptuskin 1993). This test particle calculation 
indicates, as intuitively expected, that the spectrum of 
accelerated particles is steeper than for acceleration at a 
quasi-parallel shock. As a consequence, the nonlinear ef- 
ficiency of acceleration would also have to be assumed to 
be significantly reduced. 

Fourth, the unusually high maximum CR energy pre- 
dicted by Biermann (1993) is a direct consequence of the 
assumption of a very small and energy independent CR 
diffusion coefficient, presiimably produced by the chaotic 
turbulent motions of the medium near the shock, that for- 
mally gives very fast acceleration to extremely high ener- 
gies. Leaving aside the important question concerning the 
justification of the assumption for such shock-produced 
turbulent motions also significantly ahead of the forward 
SNR shock itself (see e.g. Elfison et al. 1994, Lucek & 
Bell 2000), one has to stress that the final result should 
in any case depend upon the values of the microscopic 
CR diffusion tensor which alone allows CRs to intersect 
the shock front many times with a subsequent increase of 
their energy. As far as the macroscopic CR diffusion due 
to chaotic gas motions (Biermann 1993) is concerned, it 
can not itself provide CR acceleration. In fact, from what 
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has been said before, even our results presented below, 
based on the usual assumption about Bohm type micro- 
scopic particle diflFusion in a quasi-parallel shock and much 
less optimistic compared with those in Biermann (1993), 
must be considered as an upper limit for the CR energy 
density which can be achieved during shock acceleration. 
This leaves open the question whether the evolving SNR 
shock system can generate a strong scattering wave field 
and allows the assumption of a favorable ratio of the mi- 
croscopic parallel and perpendicular diffusion coefficients 
over a sufficiently long time, so that the maximum in- 
dividual particle energy can be increased by something 
like an order of magnitude or more for a perpendicular 
shock (Jokipii 1987; Ostrowski 1988; Ellison et al. 1995; 
Reynolds 1998). 

As mentioned above, a strong wind modifies the envi- 
ronment of massive progenitors of SN types lb and II. It 
sweeps up the ambient gas into a thin shell surrounding 
a rarefied bubble (Weaver et al. 1977; Losinskaya 1991). 
Since the typical cavity is greater than 10 pc in size and 
contains a considerable amount of matter, it can signifi- 
cantly influence the SNR evolution and CR acceleration. 
Estimates and preliminary numerical calculations indicate 
that the SN shock propagation through the progenitor's 
wind region should generate high-energy Tr'^-decay 7-ray 
emission on a detectable level (Berezinsky & Ptuskin 1988; 
Kirk et al. 1995; Berezhko & Volk 1995Q). Yet, up to now, 
all detailed numerical investigations of CR acceleration in 
SNRs have dealt with the case of a star exploding into a 
uniform homogeneous medium. 

In this paper we present a detailed extension of the 
kinetic model for CR acceleration in SNRs (Berezhko et 
al. 1994, 1995, 1996; Berezhko & Volk 1997) to the case of 
a nonuniform circumstellar medium, spherical symmetry 
still being assumed. We study the CR acceleration and 
TT^-decay 7-ray production taking into account that the 
circumstellar medium can be strongly modified by a wind 
from the progenitor star. In doing this we assume ener- 
getic particle scattering governed by the Bohm diffusion 
coefficient in the local mean magnetic field, treating the 
shock normal to be quasi-parallel to the field direction. 
We have studied two cases typical for type lb and type II 
SNe. In the first case, a star with an initial mass of 35 Mq 
at the end of its evolution explodes as SN into the cav- 
ity created by a main-sequence (MS) 0-star wind and the 
subsequent winds during the red supergiant (RSG) and 
Wolf-Rayet (WR) phases. For the case of a SN type II we 
take the example of the explosion of a star with initial 
mass 15 Mq into the cavity created by the winds emitted 
during the MS and RSG phases. In both cases we consider 



^ Due to an error in the plot routine in Fig. 1 of Berezhko 
& Volk (1995), the 7-ray fluxes in the wind SNR cases were 
plotted with an enlargement factor of 2500 which should be 
left out 



the two different interstellar number densities 0.3 and 30 
cm~'^. 

We shall not address here the question of 7-ray produc- 
tion by electron synchrotron emission, Brcmsstrahlung, 
or the inverse Compton (IC) effect on ambient low en- 
ergy photons. Especially at very high energies exceeding 
^ 10 GeV, IC 7-rays appear to dominate the emission 
from plerions like the Crab nebula (e.g. de Jager & Hard- 
ing 1992, Atoyan & Aharonian 1996), where presumably 
a relativistic wind of electron-positron pairs from the pul- 
sar is dissipated in a circumstellar termination shock deep 
inside the SNR shell. The IC effect associated with accel- 
eration of ultrarelativistic electrons inside the SNR shell, 
or at its leading shock (along with ions considered here), 
may also contribute quantitatively to the 7-ray luminosity 
of a SNR as simple estimates suggest (Mastichiadis 1996; 
Mastichiadis & de Jager 1996). Therefore in any specific 
source this leptonic contribution to the 7-ray flux needs to 
be estimated before the hadronic Tr'^-decay 7-ray emission 
can be compared with the theoretical models presented 
below. 

We briefly describe some aspects of the model in Sect. 
2 since it was in detail described in a previous paper 
(Berezhko & Volk 1997). Sect. 3 contains the results and 
Sect. 4 includes the discussion and the conclusions. 

2. Model 

2.1. CR acceleration and SNR evolution 

During the early phase of SNR evolution the hydrody- 
namical SN explosion energy Egn is kinetic energy of the 
expanding shell of ejected mass. The motion of these ejecta 
produces a strong shock wave in the background medium, 
whose size Rg increases with velocity Vg = dRg/dt. Diffu- 
sive propagation of energetic particles in the coUisionless 
scattering medium allows them to traverse the shock front 
many times. Each two subsequent shock crossings increase 
the particle energy. In plane geometry this diffusive shock 
acceleration process (Krymsky 1977; Axford et al. 1977; 
Bell 1978; Blandford & Ostriker 1978) creates a power 
law-type CR momentum spectrum. Due to their large en- 
ergy content the CRs can dynamically modify the shock 
structure. 

The description of CR acceleration by a spherical SNR 
shock wave is based on the diffusive transport equation 
for the CR distribution function f{r,p,t) (Krymsky 1964; 
Parker 1965): 

^ = 1A,2^^_^^^ + 1A(,V)^^ + Q, (1) 

dt dr dr dr dr 3 dp ' 

where Q is the source term due to injection; r, t and p 
denote the radial coordinate, the time, and particle mo- 
mentum, respectively; k is the CR diffusion coefficient; 
u = Vs ^ w. In addition, 

Wc — w for r < Rs, Wc = w + Ca for r > Rg, 
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where w is the radial mechanical velocity of the scatter- 
ing medium (i.e. thermal gas), Ca is the speed of forward 
Alfven waves generated in the upstream region by the 
anisotropy of the accelerating CRs. In the downstream 
plasma the propagation directions of the scattering waves 
are assumed to be isotropized (e.g. Drury et al. 1989). 

The thermal matter is described by the gas dynamic 
equations 



dp 1 d , o . 



0, 



dw dw d ,^ „ , 



(2) 



(3) 



dP, 



-^+^^ + 3^ir^^)P. = Ml. - l)ca^, (4) 



dt 



dr dr 



dr 



where p, 7g and Pg denote the mass density, specific heat 
ratio and the pressure of gas, respectively, and 



Pr = 



dp 



(5) 



is the CR pressure. These gas dynamic cqiiations include 
the CR backreaction via term —dPc/dr. They also de- 
scribe the gas heating due to the dissipation of Alfven 
waves in the upstream region (McKenzie & Volk 1982; 
Volk et al. 1984); it is given by the parameter a = 1 at 
r > Rs and aa = a.t r < Rs- 

We expect that the SNR shock always includes a suf- 
ficiently strong subshock which heats the gas and plays 
an important dynamical role, also in the present case of 
nonuniform background medium. The gas subshock, situ- 
ated at r = Rs, is treated as a discontinuity on which all 
hydrodynamical quantities undergo a jump. 

We assume, that the injection of some (small) fraction 
of gas particles into the acceleration process takes place 
at the subshock, that is described by the source 

Q = QsSir-Rs). 

For the sake of simplicity we restrict our considera- 
tion to protons, which are the dominant ions in the cos- 
mic plasma. At present we only have some experimen- 
tal (e.g. Lee 1982; Trattner et al. 1994) and theoretical 
(Quest 1988; Trattner & Scholer 1991; Giacalone et al. 
1993; Bennett & Elhson 1995; Malkov & Volk 1995, 1996) 
indications as to what value of the injection rate can be ex- 
pected. We use here a simple CR injection model, in which 
a small fraction rj of the incoming protons is instantly in- 
jected at the gas subshock with a speed A > 1 times the 
postshock gas sound speed Cs2 (Berezhko et al. 1990; Kang 
& Jones 1991; Berezhko et al. 1994, 1995, 1996; Berezhko 
& Volk 1997): 



Qs = 



-Sip-Pinj), Ninj ='nNi, Pinj = \mCs2, (6) 



where N = p/m is the proton number density, and m is 
the particle (proton) mass. For simplicity, we always use 
A = 2. The subscripts 1(2) refer to the point just ahead 
(behind) the subshock. 

We assume that the Bohm diffusion coefficient is 
a good approximation for strong shocks (McKenzie & 
Volk 1982), characterized by strong wave generation (Bell 
1978). This latter question has been studied very recently 
in its nonlinear consequences, by Lucek & Bell (2000) nu- 
merically, and by Bell & Lucek (2000) using an analyti- 
cal model. The conclusion is that wave generation in very 
strong shocks should not only lead to Bohm diffusion, but 
also to an amplification of the pre-shock magnetic field 
that increases the acceleration rate, as had been specu- 
lated earlier (Volk 1984). We use here the CR diffusion 
coefficient 



Kip) = PBC/3, 



(7) 



where ps is the gyroradius of a particle with momentum 
p in the magnetic field B, c is the speed of light. This co- 
efficient differs from the Bohm diffusion coefficient in the 
non-relativistic energy region, but this difference is abso- 
lutely unimportant because of the very high acceleration 
rate at p < mc. In the disturbed region we use k = KgPs/ p, 
where the subscript s corresponds to the current shock 
position r ~ Rg. An additional factor Ps/ p was assumed 
to prevent the instability of the precursor (Drury 1984; 
Berezhko 1986b). It can also be interpreted as describing 
the enhancement of magnetic turbulence in a region of 
higher gas density. 

Alfven wave dissipation (Volk ct al. 1984) as an addi- 
tional heating mechanism strongly influences the structure 
of a modified shock in the case of large sonic Mach number 
M = Vs/Cs > \/M^, Ma = Vs/ca is the Alfvenic Mach 
number, and Ca are the local sound and Alfven speeds 
correspondingly, at the shock front position r = Rs- The 
wave damping substantially restricts the growth of the 
shock compression ratio a = P2I ps at the level a w m\^^ 
which, in the absence of Alfven wave dissipation, has been 
found to reach extremely high values a w M^/"* for large 
Mach numbers (Berezhko et al. 1996; Berezhko & Ellison 
1999). 

The dynamic equations are solved under the initial 
{t = ti) conditions: 



m = 0, 

Pg = Pg0{r,tyj), 



P = P0{r, ty, 



(8) 



which neglect the background CRs and describe some am- 
bient gas distribution modified by the wind from the pro- 
genitor star emitted during a preceding time period t^. 
The time t = corresponds to the instant of SN explo- 
sion. 

The result of a core collapse supernova, many days 
after the explosion, is freely expanding gas with velocity 
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V = r/t. The density profile of the ejecta is described 
by (Jones et al. 1981; Chevaher 1982; Chevaher & Liang 
1989) 



Pej = 



where 



Ft-\v/vt)- 



V > Vt, 



(9) 



1 [3(fc - 3)M< 



15/2 



47rfc [10(fc- 5)£;,„]3/2' 



Vt 



10(fc - 5)Es 



-,1/2 



3(fc - 3)Mej 



Mf.j is the total ejected mass. For SNRs the value of the 
parameter k typically lies between 7 and 12. The pressure 
in the expanding ejecta is negligible. 

Interaction with the ambient material modifies the 
ejecta density distribution. We describe the ejecta dynam- 
ics in a simplified manner, assuming that the modified 
ejecta consist of two parts (Berezhko & Volk 1997): a thin 
shell (or piston) moving with some speed Vp and a freely 
expanding part which is described by the distribution (9). 
The piston includes the decelerated tail of the distribution 
(9) with initial velocities v > Rp/t, where Rp is the piston 
radius separating the ejecta and the swcpt-up ISM mat- 
ter. The evolution of the piston is described in the frame- 
work of a simplified thin-shell approximation, in which 
the thickness of the shell is neglected. Behind the piston 
(r < Rp), the CR distribution is assumed to be uniform. 

The high velocity tail in the distribution (9) ensures 
a large value of the SNR shock speed at an early phase 
of evolution. It increases the CR and 7-ray production 
significantly compared with the case where all the ejecta 
propagate with a single velocity (Berezhko & Volk 1997). 

In the case of a uniform ISM it was shown that the 
CR penetration through the piston plays no important 
role for SN shock evolution and overall CR production 
(Berezhko et al. 1996). In the case of an intensive wind 
from the progenitor star, a relatively large volume around 
the progenitor is occupied by a low density bubble. The 
main amount of CRs and 7-rays are produced when the 
SNR shock interacts with this bubble. Piston and shock 
sizes are comparable and much larger than the dynamic 
scale of the system during the most effective CR produc- 
tion phase. Therefore one can expect that CR penetration 
through the piston is more important in comparison with 
the case of a uniform ISM. 

The efficiency of diffusive CR penetration through the 
piston depends on the magnetic field structure, which is in- 
fluenced by the Ray leigh- Taylor instability at the contact 
discontinuity (r = Rp) between the ejecta and the swept 
up medium that is contained in the region Rp < r < Rs- 
According to Chevalier's (1982) estimate, the energy den- 
sity of the turbulent motions created by this instability is 
determined by the value of the thermal pressure Pp at the 
outer piston surface (r = Rp + 0), and may be as large 
as et = 0.2Pp. Turbulent motions in the ionized medium 
lead to the amplification of magnetic fields. We assume 
that the magnetic field grows up to the energy density 



cb = 0.5e(. This turbulent magnetic field is presumably 
distributed over a wide range of length scales. The energy 
density of the magnetic field fluctuations which resonantly 
interact with particles in the momentum range from mc 
to Pm rnay be as large as cb/ ^^(pm/iTic), which is about 
0.1 es for a typical CR cutoff momentum = lO^mc. 
Therefore we use a Bohm type diffusion coefficient in the 
piston region, corresponding to a magnetic field strength 
of B = y/SnSPp with S = 10"^. 

Radiative gas cooling is not included in our model. 
This process becomes important at the late Sedov phase 
(Dorfi 1991), when CR acceleration becomes inefScient. 

Detailed descriptions of the model and of the numerical 
methods have been given earlier (Berezhko et al. 1994, 
1995, 1996; Berezhko & Volk 1997). 

2.2. Progenitor wind bubble 

The strong wind from the massive progenitor star in- 
teracts with an ambient ISM of uniform density po = 
LAmNn, resulting to first approximation in an expanding 
spherical configuration, which is called a bubble (Weaver 
et al. 1977; Losinskaya 1991). Here m is the proton mass 
and Nb is the hydrogen number density in the back- 
ground ISM where 10% of helium by number is assumed. 
Throughout its evolution, the system consists of four dis- 
tinct zones. Starting from the center they are: (a) the hy- 
personic stellar wind (b) a region of shocked stellar wind 
(c) a shell of shocked interstellar gas, and (d) the ambient 
ISM. 

We consider here the case of a so-called modified bub- 
ble whose structure is significantly infiuenced by mass 
transport from the dense and relatively cold shell (c) into 
the hot region (b), and by thermal conduction in the op- 
posite direction, these two energy fluxes balancing each 
other to first order. During this stage the shell (c) has col- 
lapsed into a thin isobaric shell due to radiative cooling 
(Weaver et al. 1977; Kahn & Breitschwerdt 1989). 

The sizes of each zone are determined by the ISM num- 
ber density Nh, the wind speed Vw and the stellar mass- 
loss rate M (Weaver et al. 1977), and can be written in 
the form: 

/ • \ 3/10 

Ri = 4.38 



lO-^Mo/yr^ 



1/10 / ^-3/10^ ,2/5 



PC (10) 



^2000 km/s ) \ 1 cm-3 J \ 10^ yr, 
is the inner shock radius which bounds the wind region 

(a); 

\ 1/10 



R2 = 27.5 



M 



lO-SMe/yr, 



2000 km/s 



4/10 



H 



1 cm" 



-1/10 / , X 3/5 
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is the outer shock radius which separates the sheU (c) from 
the ISM. The time ty, = corresponds to the onset of a 
steady, spherically symmetric stellar wind. 

Region (a) is characterized by a negligibly small gas 
pressure Pg, a constant speed V^, the gas density 



M 



(12) 



and an essentially toroidal magnetic field (we use its value 
near the equatorial plane) 



B = B.. 



(13) 



where i?* is the radius and O the angular rotation rate of 
the star. 

The regions (b) and (c) are isobaric with a thermal 
pressure (Weaver et al. 1977) that can be written in the 
form 



Ph = 4.7 X 10" 



M 



2/5 



lO-^Ma/yr I V2000 km/s 



H 



1 cm 



3/5 



/ tw \ '^^^ dyne 
V 106 yr ) cm2 ' 



4/5 



(14) 



The gas number density of the region (b) A^b = pb/m is 
approximately uniform and given by 



iVft = 3.8 X 10" 



M 



6/35 . ,19/35 



lO-^Mg/yr j \l cm-3 

12/35 / , ^ -22/35 



2000 km/s 



V106 yr 



cm 



(15) 



The shell (c) is much denser: ^ N\j. Therefore we can 
describe the number density distribution in regions (b) 
and (c) by the combined expression 



= (TcNo ( — 



(16) 



where Nq = IANh, (Jc^o is the peak value of the shell 
number density reached at the front of the outer shock 
(r = ii2 — 0), and ctc is the outer shock compression ratio. 

The mass and heat transport between regions (b) and 
(c) are presumably due to turbulent motions in the bubble. 
We assume that the turbulent motions generate a mag- 
netic field which grows up to the equipartition value 



(17) 



2.3. Gamma-ray production 

Tr'^-decay gamma rays are produced by energetic CR pro- 
tons in inelastic collisions with gas nuclei which gener- 
ate also neutral pions that subsequently decay. The 7-ray 



emissivity of a SNR (in units of photons/s) can be written 
as (Drury et al. 1994) 

POO POO 

Q^(e^) = 1677^ / drr^ / dpp''appZ^cNgf{r,p,t), (18) 
Jo Jp~, 



where 

= 38.5 + 0.461n2(0.01876p/mc) mb 



(19) 



is the inelastic p — p cross-section (e.g. Berezinsky et al. 
1990), Z" is the so-called spectrum- weighted moment of 
the inclusive cross-section, p^ is the momentum of a CR 
particle with kinetic energy e = e^, Ng = p/m is the gas 
number density, a = 1 — dlnn/rflneis the power law in- 
dex of the integral CR energy spectrum, and n = Airp^ f 
denotes the CR differential number density. We use the ap- 
proximation for the spectrum-weighted moment as given 
by Drury et al. (1994), limiting by the value 



= TOm{0.2, 10 



1.49-2.73a-F0.53a' 



} 



(20) 



The integral 7-ray flux at the distance d from the 
source is 



3. Results 



(21) 



Detailed investigations of CR acceleration and SNR evo- 
lution in a uniform ISM have revealed important features 
of this process (Berezhko et al. 1994, 1995, 1996; Berezhko 
& Volk 1997). For a wide range of injection rates the CR 
acceleration efficiency is very high, near unity, and there- 
fore almost independent of the injection rate. Therefore 
we use here the particular value of the injection parame- 
ter ri = 10~^, which corresponds to a moderate injection 
rate for a parallel shock. This value of rj implies an in- 
jection rate which is almost one order of magnitude lower 
than that resulting from simulations of coUisionless plasma 
shocks (Quest 1988; Trattner & Scholer 1991; Giacalone 
et al. 1993), and which also corresponds to the kinetic 
MC model (e.g. Ellison et al. 1995; Baring et al. 1999) for 
parallel shocks. Our lower value of r] effectively takes into 
account the influence of the average shock obliquity for 
a quasi-spherical SNR which expands into an ISM with 
a uniform mean magnetic field: according to Ellison et 
al. (1995), and Malkov & Volk (1995,) already at an an- 
gle ~ 45° between the upstream magnetic field and the 
shock normal the injection rate is about an order of mag- 
nitude smaller than in the purely parallel shock case, if 
particle diffusion is somewhat less efficient than the Bohm 
limit. We can not exclude that even this effective injection 
rate leads to an overestimate for the global acceleration ef- 
ficiency of the SNR. 

We also restrict our consideration to a typical set of 
values of the SN parameters: hydrodynamic explosion en- 
ergy Egn = 10^^ erg, ejecta mass M^j = lOM©, k = 10. 
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The radiative eooling of the swept-up shcU of interstel- 
lar gas leads to high shell compression ratio cTc ^ 1. The 
ISM magnetic field and the CRs which are produced by 
the outer shock can nevertheless significantly restrict this 
shell compression ratio. We use here the moderate value 
CTc = 10. In fact, the final results are quite insensitive to 
the value of Gc because the SN shock becomes an ineffi- 
cient CR accelerator before reaching the region r = Ri oi 
the peak shell density. 

3.1. Type lb supernova 

As a typical example for a type lb SN we use theoretical 
results of stellar evolution with initial mass 35Mq (Garcia- 
Segura et al. 1996). According to these calculations the 
evolution consists of three stages: a MS phase with mass- 
loss rate M = 5.56 x 10~'^Mq yr~^, wind speed Vy, = 
2000 km/s, and duration At^ = 4.5 x 10^ yr; a RSG phase 
with M = IQ-'^Mq yr-i K, = 15 km/s, At^ = 2 x 10^ yr, 
and a WR-phase with M = 2.25 x IQ-^Mq yr-\ K; = 
2000 km/s, and Atyj = 2 x 10^ yr. 

According to Eq. (11) the MS wind creates a bubble of 
size i?2 = 76.8 pc. The parameters of this modified bubble 
are Nb = 5 x 10"^ cm'^, Pf, = 5.47 x 10"^^ dyne/cm^. 

The wind emitted during the RSG phase occupies a 
region of size Rf = Afu;K; = 0.3 pc. 

The fast wind from the subsequent WR star interacts 
with this dense RSG wind. After a relatively short period 
of time the WR wind breaks through the RSG wind ma- 
terial, leaving it in the form of clouds (Garcia-Segura et 
al. 1996). Taking into account that the interaction time is 
short and that the mass of the RSG wind is small com- 
pared with the bubble mass, we neglect in a zeroth approx- 
imation the influence of RSG phase on the final structure 
of the bubble. 

Using again approximately Eqs. (11) and (15), the WR 
wind inside the MS-bubblc now creates a WR bubble with 
number density Nj, = 1.02 x 10~^ cm~'^ and a shell of size 
R2 = 56.2 pc. Formally this value Nf, is twice as large as 
the number density of the MS bubble. If we disregard in 
our approximation the stellar mass lost in the WR phase 
compared to the mass of the partially swept-up MS bub- 
ble, the density of the new WR bubble can not exceed 
that of the MS bubble shell, simply for mass conserva- 
tion. It can at best be equal if we assume that the WR 
shell is completely dissipated, i.e. smeared out over the 
WR bubble. This is in fact a quite reasonable assumption 
and we shall adopt it here. Therefore we use a simplified 
structure of the bubble, whose parameters R2 = 76.8 pc, 
Nb = 5 X 10-3 cm-3, Pb = 5.47 x 10"" dyne/cm^, 
B = 7.89 /uG correspond to the MS-bubble, together with 
a hypersonic WR wind region of size Ri = 30.8 pc, corre- 
sponding to the pressure equilibrium condition pV^ = Pb. 
This approximation to the dynamics is justified to the ex- 
tent that we can neglect the overall mass lost by the star 
in comparison with the overall swept-up interstellar mass. 




Y I I < I I < I I ' I I < I I I 

1 2 

log(r/pc) 

Fig. 1. Circumstellar gas density distribution as a func- 
tion of radial distance for the case of a SNR lb in an 
ISM with hydrogen number densities Nh — 0.3 cm""^ 
(full line) and Nh — 30 cm^'^ (dashed line), and for the 
case of a SNR II with Nh = 0.3 cm-^ (dotted fine) and 
Nh = 30 cm-3 (dash-dotted line). 

The assumed circumstellar density distribution Ng is 
shown on Fig. 1 as a function of radial distance r for all 
the cases considered. 

For the WR-star we use the parameters -R* = 3 x 
10^^ cm, n = 10"^, = 50 G, which determine the 
value of magnetic field in the wind region (a). 

As we have already mentioned in the Introduction it 
is not obvious that our model (which is strictly speaking 
only valid for the case of a quasi-parallel shock) can be 
applied to the present case of a shock that is on aver- 
age almost purely perpendicular. Basically it is not clear 
whether a sufficient number of suprathermal particles can 
be injected into a quasi- perpendicular shock (Ellison et al. 
1995; Malkov & Volk 1995). On the other hand, there exist 
circumstances which should alleviate the obstacles for par- 
ticle injection and their subsequent acceleration. One may 
assume that the magnetic field is essentially disordered on 
a spatial scale Is. small compared with the main scale r. 
Over a substantial fraction of the magnetic flux tubes the 
SN shock will then be locally quasi-parallel, with particle 
injection and acceleration starting there without difficul- 
ties. If the scale Is is larger than the diffusive length of 
suprathermal particles l{pinj) = iis{Pinj)/Vs, these par- 
ticles will be at least accelerated up to the momentum 
p*, where l{p*) ~ Ib- If ^ l{Pinj)j accelerated parti- 
cles with momenta Pinj < p < p* will gain a substantial 
part of the shock energy. In this case they can amplify 
magnetic field disturbances on progressively larger scales, 
that will allow another fraction of these particles to be 
further accelerated. Ultimately, some particles might in- 
deed reach the maximum energy, which is determined by 
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the usual physical factors for a quasi-parallel shock (see 
below). However the particle spectrum is expected to be 
softer than for a quasi-parallel shock. For the suprather- 
mal particles with energy pinj ^ mVs the diffusive length 
is about / ~ 10~^i?s. Therefore the above scenario will 
take place, if the magnetic field is initially disordered on a 
scale > 10~^r, that seems not an unreasonable assump- 
tion. 

We note, that there is some experimental evidence that 
particle (at least electron) diffusive shock acceleration in 
the progenitor wind actually takes place (see e.g. Chevalier 
1982, Kirk et al. 1995, and references there). 

We shall use here the assumption of Bohm type CR 
diffusion for a locally quasi-parallel shock. It implies an 
acceleration rate that is almost independent of the mag- 
netic field structure. However, in the light of our above 
discussion of this field structure we consider our results as 
upper limits to the real CR and gamma-ray production 
efficiency. 

Starting from some initial instant ti = Rpi/Vpi = 
6.34 yr after the SN explosion which corresponds to an 
initial piston position Rpi — 3 x 10^^ cm and a piston 
speed Vpi = 15 x 10^ km/s, we explore the solution of the 
nonlinear dynamic equations following the SN shock prop- 
agation through the successive bubble zones. Note that 
the value of the initial position Rpi from which we start 
our calculations has as little physical importance as the 
short subsequent period of SNR evolution during which 
the quasi-stationary character of system evolution is es- 
tablished. 

Fig. 2 illustrates the CR characteristics, SN shock pa- 
rameters and expected 7-ray spectra produced during SN 
shock propagation through the supersonic WR wind re- 
gion r < i?i . 

In Fig. 2a the calculated CR distribution function at 
the subshock fs{p,t) — f{r — Rs,p,t) is presented for 
five different instants of time. The shape of the CR spec- 
trum is determined by the fact that the shock is essen- 
tially modified by CR backreaction. Very soon after the 
beginning of CR acceleration the total shock compression 
ratio a = pij ps reaches about a = 16.5 and then slowly 
decreases with time. Although the total compression ra- 
tio is much higher than predicted by the hydrodynamical 
two-fluid model (Jones & Kang 1992), the subshock never 
disappears and its compression ratio remains nearly con- 
stant at the level ~ P2/P1 ~ 3.4 (sec Fig. 2b). At small 
momenta Pinj < p^mc the CR distribution function is 
almost a pure power-law fs oc with the power-law 
index 

qi = 3a'J{ai - 1), (22) 

determined by the effective subshock compression ratio 
a'g = as{l — 1/Moi). At momenta 10 mc;^p;^Pm, where 
Pm ^ lO^mc, the distribution function is also close to a 
power-law form fg oc with the index value q = 3.6. 



It is important to note, that our results confirm the 
conclusion that the accelerated particle spectrum can not 
be harder than / oc p~^'^ even for a very large shock com- 
pression ratio (7 (Berezhko 1996; Malkov 1997; Berezhko 
& Ellison 1999a), whereas the test particle approach (see 
expression (27) below) in the case ct » 1 predicts the 
limiting value of the power law index q = 3. 

An interesting point is that in the wind the CR dis- 
tribution functions fs{p) at different evolutionary phases 
are almost selfsimilar to each other except during the ini- 
tial short period i < 10 yr when the system undergoes 
the transition from the initial conditions at t = ti to 
the quasi-stationary state. It is a result of the fact, that 
during the propagation through the low density super- 
sonic wind the shock speed remains nearly constant. The 
analysis of the dynamical equations shows that in this 
case the CR distribution function has a selfsimilar form 
f{r,p,t) = ^{r/Rs,p)/t^ (see Appendix) which is roughly 
consistent with the numerical results presented on Fig. 2a. 
It also explains the time constancy of the maximum par- 
ticle momentum obtained by Volk & Biermann (1988), as 
mentioned in the Introduction.. 

A more detailed consideration includes the shock de- 
celeration. In the case of the ejecta density described by 
the selfsimilarity distribution Eq. (9) at the early phase, 
when Vs > ft, the expected expansion law is 

(23) 

with ly — {k — 3)/{k — 2) (Chevalier 1982), which in our 
case has the value v = 0.875. The decrease of the shock 
velocity 

Vs=u^<x t"-^ (24) 

leads to a decrease of the injection momentum pi^j ~ 
mVs oc t''~^ that introduces an additional time depen- 
dent factor in the CR distribution function / cx pj^J'^ ^ 
i(!^-i)(9<.-3) a,t least at nonrelativistic energies. It results 
in the dependence / oc t~^^ where p, = 2u — {u — l){qs — 3) . 
In our case k = 10, as = 3.4 this gives qs = 4.25 and 
p = 1.91. Due to the selfregulating property of the non- 
linear acceleration process, CR particles absorb some con- 
stant fraction of the shock energy; therefore at relativistic 
energies oc PsV^ oc t~'^. 

A direct consequence of this fact is the evolution of the 
overall, i.e. volume integrated CR spectrum 

/•OO 

7V(p,<) = 167rV / drr^f{r,p,t), 
Jo 

which is shown in Fig. 2c as a function of momentum 
p for five subsequent evolutionary phases. The shape of 
the spectrum and the value of the cutoff momentum pm 
change slowly during the evolution. At high energies 
10^<p/mc<10^ the overall CR spectrum is close to the 
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Fig. 2. CR and gas dynamics in a WR wind region. CR distribution function at the shock front (a), volume integrated 
overall CR spectrum (c), as functions of momentum, and overall integral 7-ray spectrum normalized to a distance of 

1 kpc (c), as function of 7-ray energy, at five different times t. Total shock compression ratio cr, subshock compression 
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form iV cx p~^-^ . It is interesting to note that the overall 
CR spectrum is even harder than the local CR spectrum 
at the shock front Hg = Airp^fsip), because the size of the 
upstream region / = Ks{p)/Vs occupied by CRs with mo- 
mentum p is proportional to p. At relatively low momenta 
l{p) « Rg — Rp and the overall spectrum N{p) mainly 
consists of the CRs situated in the downstream region 
Rp < r < Rs] therefore in this case N{p) oc ns{p). At high 
momenta p ^ Pm diffusive length l{p) becomes compara- 
ble with the downstream region size Rg — Rp ~ _Rs/(3fT)' 
therefore the overall spectrum N{p) becomes progressively 
harder with increasing p due to upstream region. The 
amplitude of the spectrum grows with time because the 
volume occupied by accelerated CRs increases with time 
V (xR^g that gives N oc Vf^p'^ oc. t^"-!". The index Su - /j, 
varies from 0.72 at p < mc to 0.63 at p > mc in good 
agreement with the numerical results. 

One can compare the calculated and the expected 
value of the CR cutoff momentum, which is determined 
by the expression (Berezhko 1996) 



Pm RsiVg - Vn, 



mc AKs{mc) 
where 



A = [2 + 2b + e + d - {i^ - l)/p]qt/{5 - Qt), 



(25) 



(26) 



and the parameter d = (r'Vw)2<j/[{o- — 1)14] describes 
the effect of particle adiabatic cooling in the downstream 
region, and the dimensionless parameters 

v = din Rg/ dint, b = din fg/dlnRg, e = vdlnKg/dhiRg 

describe the time variation of the shock radius, of the CR 
distribution function and of the CR diffusion coefficient, 
respectively. In addition 



qt = 3a' /{a' - 1) 



(27) 



is the lower limit for the power law index, and a' = 

(Vg - K, - Ca)/u2 = cr(l - 1/Ma) is thc cffcctivc total 
shock compression ratio which includes the effect pro- 
duced by the outward propagation of Alfven waves in 
the upstream region r > Rg. By definition thc value of 
Pm is that momentum where the local power law index 
q = —d{lnf)/d{lnp) drops to the value the value q = 5. 

In thc case under consideration thc gas velocity w is al- 
most constant in the downstream region, which gives d = 2 
(Berezhko 1996). Taking into account that Rg oc 
fs cx Kg OC Rg wc have u = 0.875, b = —2, and e = 
0.875, and the expression for the expected cutoff momen- 
tum can be written in the form 



Pn 

mc 



= 2.5 X 10"* 



R, 



50 G 



3 X 1012 cm 



2 X 10^ km/s 



10 yr 



(28) 



-0.125 



which is in a good agreement with thc numerical re- 
sults (see Fig. 2a,c). Note that the main factors which 
determine the value of Pm are the adiabatic CR cooling 
in thc expanding medium, thc finite time increase of the 
shock size, and the decrease of the shock velocity, but not 
the time factor, as it is frequently assumed (e.g. Lagage 
& Cesarsky 1983). 

During the time interval t < 5922 yr of shock prop- 
agation through region (a) of the supersonic wind CRs 
absorb only a small fraction of the explosion energy E,, ~ 
0.008i?s„ (see Fig. 2f) due to the small amount of swept 
up matter Mgw = O.SMq compared to the ejected mass 
Mej = WMq, where 



Mg. 



4-77 



drr /9, 



(29) 



p{r) is the wind density, determined by the expression 
(12). 

The 7-ray spectrum is shown in Fig. 2c. It is ex- 
tremely hard in the energy range 10^° ^ ;$ lO^^ eV dur- 
ing the early phase t < 500 yr and becomes progressively 
steeper as time proceeds. This bcihavior is a consequence 
of the fact that the 7-ray flux F^ = F!.^ + i^" consists of 
two different parts. The first, F^, is due to CR interaction 
with the swept-up matter which lies between the piston 
and shock surfaces and has the mass Mg^. In the case we 
consider here, the CRs are almost uniformly distributed in 
thc downstream region Rp < r < Rg. Therefore we have 
approximately (Berezhko & Volk 1997) 



F'. oc M,. 



,ec2P: 



(30) 



where ec2 is the CR energy density in the downstream 
region. In our case €02 oc PgV^, Mg^ = MRg/Vw, 7 ~ 1.7, 
giving thc expected dependence F!^ oc t~^'^ which is in a 
agreement with the numerical results at f > 10 yr. The 
calculated 7-ray flux can be represented in the form 



F'{1 TeV) 



Flo 



d 



1 kpc 



t 



10 yr 



-1.1 



(31) 



with FIq 



2 X 10~i2 cm~2 s~i, which is in rough agree- 
ment with the estimates of Bcrczinskii & Ptuskin (1989), 
and Kirk et al. (1995) which do not include CR pene- 
tration into the ejecta with subsequent 7-ray generation. 
The energy spectrum of this component F^ oc e;:^'^'^ is 
steeper than the CR integral spectrum, because upstream 
CRs interacting with the relatively low density medium 
contribute much less to the 7-ray spectrum than to the 
overall CR spectrum 

The second part of 7-ray flux i^" has its origin in thc 
ejecta material whose CR energy density is Ccs- Therefore 
we have 



cx Me,e,3p2 
The ratio of thc two components 
F^'/F; = {M,^/Mg^){e,3/e,2) 



(32) 
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Fig. 3. Integral flux of 7-rays with energy e-^ >1 TeV, 
normalized to a distance of 1 kpc, as a function of time 
for the cases of a SNR lb with = 0.3 cm'^ (full 
line), SNR lb with Nh = 30 cm'^ (dashed Une), SNR 
II with = 0.3 cm-3 (dash-dotted line), and SNR II 
with Nil = 30 cm" 3 (dotted line). 



is determined by the two factors M^jjMsw and ecil^c2- 
The swept up mass M^^ = MRs/V^ increases with time 
and reaches the value = 0.03M<y at t ^ 6000 yr. At 
t < 100 yr the ratio 603/602 grows from zero to about 0.05 
and then remains nearly constant. Therefore, during the 
initial period t < 10^ yr, the 7-ray flux at e-^, ^ 1 TeV 
is dominated by the second component whereas at t > 
6000 yr the first component becomes dominant at almost 
all energies. 



t 



Since 6c3 oc psV^ oc t ^ 
and numerically 



F''{1 TeV) = 4.6 x 10 



d 



-12 



1 kpc 



a.t t > 100 yr we have F" ex 



(33) 



100 yr 



The second component F"{e^) is much harder than 
F^{ej) (see Fig. 2e) because CRs with higher energy pen- 
etrate into the piston more effectively. Therefore the ex- 
pected flux F^{e^) becomes progressively steeper as the 
contribution of F" decreases with time (see Fig. 2e). 

The time dependence of the expected integral flux of 
TeV 7-rays at the distance d = 1 kpc is shown in Fig. 3. 
It also includes the epochs when the SNR shock leaves the 
wind region and enters the bubble and the shell regions, 
to which we turn now. 

For t > ti = 5922 yr the SN shock propagates through 
the hot region (b) of the MS bubble. When the SNR shock 
intersects the boundary r = R\, which is also a strong 
discontinuity, a new secondary shock arises (Shigeyama & 
Nomoto 1990). It propagates several times between the 



SNR shock front and the piston surface, and compresses 
and heats the medium. It thus provides the transition to 
the new quasi-stationary state which corresponds to the 
SNR shock propagating through region (b). We neglect 
this complicated transition phenomenon. We rather de- 
scribe the beginning of the SN shock propagation in the 
region (b) in the following simplified manner. We start 
with a pure gas shock of size Rs = Rp + 0, piston size 
Rp = J?i and piston speed Vp = Vpi = 4127 km/s which 
is the final piston speed in the region (a). We neglect the 
CRs produced during the previous period t < ti. This un- 
derestimates the CR and 7-ray production in region (b). 
But the influence of previously produced CRs is not very 
large because the CR production rate increases sharply 
when the SNR shock intersects the boundary r = on 
account of the much higher gas density in region (b) (see 
Fig.l). 

Fig. 4 illustrates the dynamics of the system during 
this period. 

Note that the CR distribution functions in Figs. 2a and 
4a are measured in the same units; this is also the case for 
the overall CR spectra in Figs. 2c and 4c. One can see that 
already at t = 6100 yr the CR number density and the CR 
distribution function increase by more than two orders of 
magnitude compared to the time t = t\ = 5922 yr. In the 
region (b) efficient CR acceleration lasts less than four 
thousands years: only for t^lO^ yr the shock produces 
a CR spectrum that is relatively hard in the relativistic 

:comes pro- 



it be 



energy range (see Fig. 4a). For t > 10"* 
gressively steeper so that the CR number density at energy 
e = 10 TeV, which generates TeV-energy 7-rays, has al- 
ready at t = 3 X 10^ yr decreased by more than two orders 
of magnitude. Due to this reason TeV-energy 7-rays at 
t > 10* yr are produced by CRs accelerated at the earlier 
stage t^W^ yr, and the contribution of freshly acceler- 
ated CRs to the 7-ray production is negligible. The 7-ray 
flux decreases with time for t ^ 10* yr because of adiabatic 
cooling of the CRs accelerated during the previous stage. 

Due to its high temperature, region (b) is charac- 
terized by a relatively large value of the sound speed 
Cs = 188 km/s. When it breaks into the region (b) the 
shock has a velocity of about Vs = 4500 km/s. Therefore 
even the initial Mach number M = 24 is not very high. 
Only during an initial period 10* yr the shock compres- 
sion ratio (T is substantially larger than 4 (see Fig. 4b), 
and high energy CRs are produced with relatively high 
efliciency. At later stages 10* yr the shock decelerates 
rapidly (Fig. 4d) and becomes too weak to accelerate high 
energy CRs efficiently. Due to the small Mach number the 
total shock compression ratio a at t^lO* yr is less then 
4 and even the lower limit qt for the power law index at 

10* yr becomes greater than 5 (Fig. 4b). The CR spec- 
trum fs oc produced by the SN shock at this stage 
is very steep {q > qt > 5) which is tantamount to a low 
acceleration efficiency. 




Fig. 4. The same as Fig. 2 but for SNR sliock propagation through the WR bubble which was reached after 5922 yr. 
The normalization parameters in Fig. 4c are, respectively, Ri = 30.8 pc, Vpi = 4127 km/s. 
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One can see from Fig. 4c that the number of high en- 
ergy CRs withp/mc > 10^, produced in the SNR, increases 
only during the initial period 10* yr and then remains 
nearly constant during shock propagation in the bubble. 
The lower energy particles, on the other hand, continue 
to be accelerated beyond this time, steepening the over- 
all spectrum at lower energies, whereas Ec remains nearly 
constant. At the end of this stage about 22 percent of the 
explosion energy is transformed into CRs (see Fig. 4f). 
Since the chemical composition of the bubble gas is more 
and more modified towards a pure interstellar composition 
the further out in radius the shock extends, we conclude 
from Fig. 4c that the chemical composition of the CR spec- 
trum becomes more wind material-like (heavier) towards 
higher particle energies. The same tendency is seen in the 
later SN II case (see Fig. 5c). This effect is peculiar to wind 
SNe but it is only one of several factors which lead to an 
increasingly heavier chemical composition of the Galactic 
CRs with increasing particle energy. 

Approximately, neglecting the mass in the wind, the 
behavior of CR and 7-ray production in the bubble is de- 
termined by the dynamical scale length Iq and the corre- 
sponding time scale 



to = Iq/Vo, where Vo = ^^2Esn/Mej 

is the mean cjccta velocity. As in the case of a uniform 
medium Iq is the length over which the amount of swept- 
up material equals the ejected mass M^j. The principal 
difference to the case of a uniform medium is that in the 
case of the bubble the shock expansion starts from the 
initial size Rg w Ri which is much larger than Iq. 

Therefore the appropriate definition for the dynamic 
scale in the case under consideration is 



lo 



M. 



inRipb' 



(34) 



where pb is the bubble density. It gives the value = 
4.92 pc which is much smaller than Ri. The corresponding 
time scale is to = 1522 yr. One can see from Fig. 3 that 
in agreement with this time scale the TeV-encrgy 7-ray 
emission reaches its peak value 



10-13 cm-^s-i 



at 



t Ri fi -|- 2fo ~ 10* yr, that is at the beginning of the Sedov 

phase as in the case of a uniform ISM. Later on the 7-ray 
flux gradually decreases (see Fig. 3) due to the decrease of 
the CRs production because of the decrease of the shock 
strength (sec Fig. 4b). 

For a long time (at least up to f = 8 x 10^ yr) the CR 
energy content Ec remains nearly constant. In the case 
under consideration the adiabatic cooling of CRs is less 
important compared with the uniform ISM case because 
the shock size at the Sedov phase varies over a very small 
range: during the period from i « 7 x 10^ yr to i w 7 x 
10** yr the shock size increases only by a factor of two (see 
Fig. 4d). 

At i > 10* yr the CR energy content at the shock 
front is mainly in the form of low energy particles because 



freshly accelerated CRs arc characterized by a progres- 
sively steeper spectrum (see Fig. 4a). Note that the local 
peak (bump) in p*/s(p) at p ~ lO^mc represents the so- 
called escaping particles (Berezhko 1986a; Berezhko et al. 
1996). Fig. 4c shows that the overall CR spectrum N{p, t) 
in the relativistic energy range remains nearly constant 
for t > 10* yr. It has an almost pure power law form 
(X p~^^^ in the momentum interval lO^p/mc^lO^ 
(Fig. 4c). 

The efficiency of high energy 7-ray production in- 
creases again for i > 2 x 10* yr due to the increasing 
gas density (see Fig. 4e). But even at f ~ 10^ yr, when 
Ng{Rs) ^ 1 cm~^ the TeV-energy 7-ray flux remains sig- 
nificantly lower than 10"^^ cm-^s"*. 

Note that during this stage the total number of high 
energy CRs remains almost constant since the shock be- 
comes too slow and weak, and does not accelerate ener- 
getic particles any more. Nevertheless one can expect a sig- 
nificant increase of high energy 7-ray production later on 
when the strong population of previously produced CRs 
will expand outward and finally reach the shell (c) and its 
boundary. During this stage an almost constant number 
of CRs progressively interacts with an increasing amount 
of gas and this will lead to an increasing 7-ray production. 
Unfortunately, it is not simple to describe the dynamics of 
the system during this late stage. One has for example to 
expect a strong increase of the CR diffusion coefficient due 
to the decrease of the CR energy density as well as magne- 
tohydrodynamic wave damping by ion neutral friction in 
the partially ionized shell gas. This suggests a much faster 
outward loss of CRs compared to the Bohm diffusion case, 
unless the SNR is surrounded by a hot, fully ionized gas 
which reflects the particle back across the shell. 

In Fig. 3 we also present the time dependence of the 
TeV-energy 7-ray flux for the same case of a SNR lb 
but for a higher (asymptotic) ISM number density Nh = 
30 cm~^. In this case we have Ri = 14.6 pc, R2 = 34.7 pc, 
Nb = 2.82 x 10-2 cm-3, P(, = 1.1 x 10"" dyne/cm^, 
Bb = 16.7 fxG. One can see that the character of the 
SNR evolution is very similar to the previous case. For 
t < 2544 yr SN the shock propagates through the region 
(a) of the supersonic WR-wind and its evolution is initially 
identical to the previous case. 

The bubble region, reached by the SNR shock at 
ti = 2544 yr, is characterized by the length scale lo = 
1.8 pc and the time scale to = 555 yr. The value = 
3 X 10-13 cm-^s-i, reached at t h + 2to = 3654 yr, 
would be the peak flux for a uniform medium with num- 
ber density A^o = Nb- In contrast to the previous case, the 
7-ray flux continues to increase for t > 3654 yr because at 
this period the SN shock enters the region (c), where the 
gas density starts to increase rapidly (see Fig. 4d). 

One can see from Fig. 3 that the 7-ray flux reaches the 



value F-v 



10-^° cm-^s-i at i Ri 3 X 10* yr. A compari- 
son with the previous case shows that this peak flux scales 
as Njf as in the case of a uniform ISM (e.g. Berezhko & 
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Volk 1997). However, it is essentially lower (at least by a 
factor of 100) than in the case where the ISM is uniform 
with the same density Nh- The main reason of this low 
efficiency of 7-ray production is that in the case under dis- 
cussion the majority of CRs is produced in the low density 
bubble. When the SN shock enters the high density shell 
region (c) it becomes weak and therefore an inefficient CR 
accelerator. In this phase, in comparison with the case of a 
uniform ISM, the same amount of CRs mainly produced 
at previous evolutionary phases generates a much lower 
7-ray flux because of the poor spatial overlap between the 
CR and gas distributions (see Berezhko & Volk 1997). 
The peak of the gas density distribution lies just behind 
the shock front, whereas the CRs are sitting deeper inside 
where the density of the gas is much lower. 

3.2. Type II supernova 

We model the type II SN case as a progenitor star with ini- 
tial mass 15A'/q (e.g. Leithcrcr ct al. 1992) that spends a 
time period At^ = 4 x 10^' yr on the main-sequence with a 
mass-loss rate M = 2.5 x 1O^^M0 yr"-'^ and wind velocity 
Vw = 2000 km/s, and then the time Ai^ = 10^ yr in the 
RSG phase with mass-loss rate M = 2 x 1O~^M0 yr~^ and 
wind velocity 14; = 15 km/s. The MS wind creates a bub- 
ble of size i?2 = 61 pc in the ISM with Nr = 0.3 cm"^. Ac- 
cording to Eqs. (14) and (15) the bubble is characterized 
by AT;, = 6.6 X IQ-^ cm"^ and = 4.4 x lO'^^ dyne/cm^. 

The RSG wind occupies the region of size Rf = 
Vyjl^tw = 1-54 pc. At this point the ram pressure pV^ 
of RSG wind exceeds the thermal pressure Pb in the bub- 
ble. Therefore we neglect the shell which can be formed 
in the RSG wind due to its interaction with the ambient 
bubble material. We model the transition zone between 
the RSG wind and the bubble by the set of parameters 
p = p{Rf){Rf/rf-^, = V^{Rf){R}/rf. We introduce 
this zone to match smoothly the gas densities Ng between 
regions (a) and (b), sec Fig. 1. It contains a small amount 
of gas and plays no role in the overall SNR evolution. 

We use a magnetic field strength _B = 2 x 10~* G in 
the RSG wind at the distance r = 10^^ cm. It formally 
correspond to = 1 G, i?* = 3 x 10^^ cm and Q = 
3 X 10~®. Our B-field is about 60 times smaller than that 
assumed by Volk & Biermann (1988). Their magnetic field 
implies an Alfven velocity Ca which is considerably larger 
than the wind speed Vw This renders the existence of 
the RSG mass outflow as a supersonic wind problematic 
(Axford 1994). In our case the RSG wind is superalfvenic 
with Ca w 0.5Vu,. It might be possible to model a slow 
outflow starting with a much higher stellar magnetic field 
which is then subalfvenic. However, we shall not attempt 
such a dynamical construction in this paper. 

We start the computation of the SNR evolution with 
the initial conditions ti = 3.17 yr, Rpi = 10^^ cm, Vpi = 
10* km/s. 



Fig. 5 illustrates SN shock evolution, CR and 7-ray 
production for an ISM with Nh = 30 cm""^. The charac- 
ter of the CR spectrum fs{p,t), shown in Fig. 5a, and of 
the overall CR spectrum N{p,t), plotted in Fig. 5c, are 
similar to the previous case of the SN lb for t < 500 yr 
when the SN shock propagates through the region (a) of 
the supersonic wind. The value of the cutoff momentum 
Pm ~ 2 x 10^ mc is consistent with formula (28). The ac- 
celeration efficiency and the corresponding shock modifi- 
cation arc very high. During the initial stage i « 10 yr the 
shock compression ratio reaches the value cr = 15 and then 
slowly decreases due to the shock deceleration (see Fig. 
5b) . As a consequence the high energy part of the CR spec- 
trum is extremely hard N rx p~^'^ at 10'' <p/toc < 10^. 

It is interesting to note that the 7-ray spectrum oc 
^-0.5 (Yig^ is niuch steeper at 10"'^"'^ < e-, < lO"'^"' cV than 
the integral CR spectrum: the CRs in the upstream region 
r > Rs, which make the integral CR spectrum extremely 
hard, play a much less important role in the formation 
of the 7-ray spectrum, because they occupy a relatively 
lower density region. 

In the present case the 7-ray generation within the 
ejecta is negligible. Therefore, according to relation (30), 
the expected 7-ray flux F^(l TeV) should be about 5 x 10^ 
times larger than the flux F! (l TcV) in the case of the 
WR wind. One can see from Fig. 3 and Fig.5e, that the 
calculated 7-ray flux can be represented in the form (31) 
with Fio = 1.2 X 10~^cm~^s~^ in agreement with relation 
(30). 

During these initial 500 years the SN shock sweeps 
up Msw = '^Mq, CRs absorb about 10% of the explosion 
energy (see Fig. 5f) and the expected 7-ray flux -^^(l TeV) 
exceeds the value 10"^" cm~^s~^ (see Fig. 3). 

During the period 500 < t < 1000 yr the SN shock 
propagates through the intermediate zone (see Fig. 5d) 
which contains only a small amount of matter. The CR 
energy content slowly increases, but the CR spectrum be- 
comes progressively steeper. Together with the decrease of 
the gas density this leads to a decrease of the 7-ray flux 
F^ a (see Fig.3). 

For t > 10^ yr the SN shock propagates through the 

bubble whose characteristics are similar to the bubble 
around the progenitor of a type lb SN, in an ISM with 
the same density Nh (see Fig. 1). Therefore at this stage 
the expected 7-ray production continues to decrease up 
to the time t = 4 x lO"' yr, when the 7-ray production by 
CRs accelerated during previous stages drops to a level 
that corresponds to a flux F^. ^ lO"""^^ cm~^s~-'^ typical 
for a uniform medium with A^o = A^fa- For t > 4 x 10"^ yr 
7-rays produced in the bubble material start to dominate 
and the 7-ray flux increases with time. At t = yr 
the SN shock reaches the shell region (c) where the den- 
sity increases. Therefore the expected 7-ray flux increases 
more rapidly for t > 10^ yr (see Fig. 3), even though the 
CR production becomes quite low in this phase (see Fig. 
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5f ) . More than 40% of the explosion energy goes into CRs 
during the full SNR evolution. 

In Fig. 3 we present also the calculated 7-ray flux ^^(1 
TeV) for the case Nh = 0.3 cm~^. The time-dependence 
of i^-y(l TeV) for < < 2 X 10^ yr is identical to the previous 
case. For t > 2x 10^ yr the SN shock propagates through 
the bubble whose characteristics are now similar to the 
case of a type lb SNR with Nh = 0.3 cm"-^. Therefore, the 
expected 7-ray production during late phases t > 10* yr 
is close to that case (see Fig. 3). 

4. Summary 

Our numerical results show that when a SN explodes into 
a circumstellar medium strongly modified by a wind from 
a massive progenitor star, then CRs are accelerated in the 
SNR almost as effectively as in the case of a uniform ISM 
(Berezhko et al. 1994, 1995, 1996; Berezhko & Volk 1997): 
about 20 40% of the SN explosion energy is transformed 
into CRs during the active SNR evolution. 

During SN shock propagation in the supersonic wind 
region very soon the acceleration process reaches a quasi- 
stationary level which is characterized by a high effi- 
ciency and a correspondingly large shock modification. 
Despite the fact that the shock modification is much 
stronger than predicted by a two-fluid hydrodynamical 
model (Jones & Kang 1992), the shock never becomes 
completely smoothed by CR backreaction: a relatively 
strong subshock always exists and plays an important dy- 
namical role. As in the case of a uniform ISM, the spurious 
complete shock smoothing in hydrodynamic models is the 
result of an underestimate of the role of geometric factors. 

Due to the relatively small mass contained in the su- 
personic wind region CRs absorb there only a small frac- 
tion of the explosion energy (about 1% in the case of a 
SN type lb, and 10% in the case of a SN type II) and 
the SNR is still very far from the Sedov phase after hav- 
ing swept up this region. Therefore we conclude, that the 
CRs produced in this region should not play a very signif- 
icant role for the formation of the observed Galactic CR 
energy spectrum. 

The peak value of the CR energy content in the SNR 
is reached when the SN shock sweeps up an amount of 
mass roughly equal to several times the ejected mass. This 
takes place during the SN shock propagation in the mod- 
ified bubble. In a purely adiabatic bubble the sweep-up 
would occur at the beginning of the shell. Compared with 
the uniform ISM case the subsequent adiabatic CR de- 
celeration is less important in the case of a modified cir- 
cumstellar medium. The main amount of CRs in this case 
is produced when the SN shock propagates through the 
bubble. In this stage the dynamical scale length is much 
smaller than the shock size. Therefore the relative increase 
of the shock radius during the late evolution stage and the 
corresponding adiabatic effects are small. These configu- 
rational properties lead to potentially interesting changes 



of the CR chemical composition with particle energy (see 
Sect. 3.1) 

In the case of the modified circumstellar medium 

the CR and 7-ray spectra are more variable during the 
SN shock evolution than in the case of a uniform ISM 
(Berezhko & Volk 1997). At the same time the form of 
the resulting overall CR spectrum is rather insensitive to 
the parameters of the ISM as in the case of uniform ISM. 
The reason is that the main amount of CRs are produced 
in the latest phase which has still a strong enough shock. 
Roughly speaking, the overall CR spectrum (except the 
most energetic CRs) is mainly formed at the stage when 
the shock compression ratio lies between 4 and 5. 

The maximum energy of the accelerated CRs reached 
during the SNR evolution is about 10^* eV for protons in 
all the cases considered. 

Our results confirm the important conclusion, reached 
for the case of a uniform ISM before, that the diffusive 
acceleration of CRs in SNRs is able to generate the ob- 
served CR spectrum up to an energy ~ lO^** eV, if the 
CR diffusion coefficient is as small as Bohm limit. This 
disregards the possibility of turbulent field amplification, 
as discussed in Sect. 2. 

In the case of a SN lb the expected TeV-energy 7-ray 
flux, normalized to a distance of 1 kpc, remains lower than 
10"^^ cm~^s"^ during the entire SNR evolution if the ISM 
number density is less than 1 cm~^ except for an initial 
short period t < 100 yr when it is about 10""'^^ cm~^s~^. 
Only for a relatively dense ISM with Nh = 30 cm~"^ the 
expected 7-ray flux is about 10"^*^ cm~^s~^ at late phases 
t > 10"* yr. A similar situation exists at late phases of SNR 
evolution in the case of SN II. It is interesting to note that 
the expected 7-ray flux is considerably lower, at least by a 
factor of himdred, compared with the case of uniform ISM 
of the same density Nh- This confirms the preliminary 
result reported earlier (Berezhko & Volk 1995) 

The type II SN explodes into the dense wind of the 
red supergiant progenitor star. During the first several 
hundred years tm after the explosion, the expected TeV- 
energy 7-ray flux at a distance d = 1 kpc exceeds the 
value 10"^ cm~^s~^ and can be detected up to the dis- 
tance dm = 30 kpc with present instruments like HEGRA, 
Whipple or CAT. This distance is of the order of the di- 
ameter of the Galactic disk (see also Kirk et al. 1995). 
Therefore all Galactic SNRs of this type whose number is 
Nsn = Vsntm should bc visiblc. But in this case we can 
expect at best N^n ~ 10 such 7-ray sources at any given 
time. 

The typical value of the cutoff energy of the expected 
7-ray flux is about 10^^ eV, if the CR diffusion coefficient 
is as small as the Bohm limit. In this respect the negative 
result of high-threshold arrays (Borione et al. 1995; Allen, 
G.E. et al. 1995; Allen, N.H. et al. 1995) in searching of 
7-ray emission from Galactic SNRs is not surprising be- 
cause their threshold Eth ^ 50 TeV exceeds the cutoff en- 
ergy of the expected 7-ray flux; marginally this also holds 
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for the negative results of the lower threshold Eth > 20 
TeV AIROBICC array (Prahl & Prosch 1997; Prosch et 
al. 1996). It is less obvious how to interpret the nega- 
tive results of imaging atmospheric Cherenkov telescopes 
with thresholds less than about 1 TeV (Mori et al. 1995; 
Lessard et al. 1997; Hess et al. 1997; Hess 1998; Buckley 
et al. 1998). For core collapse SN of types II or lb with 
quite massive progenitors one can in part explain this fact 
by the extremely low Tr'^-decay 7-ray intensity expected 
from such SNRs during the period of SN shock propaga- 
tion through the low-density hot bubble. An alternative 
possibility relates to the assumption of the Bohm limit for 
the CR diffusion coefficient which can be too optimistic, 
in particular for the quai-perpendicular geometry in wind- 
blown regions from rotating stars. For a slightly more gen- 
eral discussion of SNR 7-rays in stellar wind cavities, see 
Volk (1997). 
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Appendix A: Similarity solution 

Consider a shock of radius Rg = Vgt that expands with 
constant speed Vg into the wind region, whose parameters 
are described by expressions (12), (13). We introduce the 
similarity variables 

X = r/Rs, 

f{r,p,t)=^{x,p)/e, 
P,{r,t)=IL,{x)/t'', 
w{r,t) = W{x), 

p{r,f) = n{x)/t\ 
Pg{r,t)=Ilg{x)lt\ 

Then the diffusive transport equation for the CR distri- 
bution function, eq.(l), and the gas dynamic equations 
(2)-(4) can be written in the form 



{W-xV,) 



dx 



2|l:,.^»')n, = <..(..-i)c.^, 



2^-\-x 



x^K- 



dx x'^V, dx dx V, dx 



-+ 



1 ^x^H^)f^ + !fc|^^(.-l), 



xWg dx 3 dp ATrmpf^jVs 

^{{xVs - W)Q] - - 3V,n = 0, 

dW d 

ixV.-w)n—^-iii, + u.), 



taking into account that in the wind region the Alfven 
speed Ca is constant and small compared with the shock 
speed Vs, and plausibly assume that the CR diffusion coef- 
ficient has the form n{r,p,t) = K{x,p)Rs- The boundary 
conditions do not contain the time explicitly. Therefore 
the above similarity solution is appropriate. The only fac- 
tor which violates these assumptions it is the initial con- 
dition for CRs, which contains the time t = U in explicit 
form. Therefore the exact solution will deviate from the 
similarity solution only during some short initial period of 
several ti, as long as the shock speed is constant. 
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Fig. 5. The same as Fig. 2 but for the case of the remnant of a type II SN in an ISM with the hydrogen number 
density Nh = 30 cm~'^, and encompassing the entire SNR evolution. The normahzation parameters are, respectively, 
given by Rpi = 10^^ cm, Vpt = 10^ km/s 



